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Abstract. Precision radial velocity data for HD 208487 has been re-analyzed using a new Bayesian 
multi-planet Kepler periodogram. The periodgram employs a parallel tempering Markov chain 
Monte Carlo algorithm with a novel statistical control system. We confirm the previously reported 
orbit (Tinney et al. 2005) of 130 days. In addition, we conclude there is strong evidence for a second 
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planet with a period of 998+^ davs ' an eccentricity of 0.19+q % and an Msin/ = 0.46+qJ|M/. 
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High-precision radial-velocity measurements of the reflex radial velocity of host stars 
have lead to the discovery of over 150 planetary candidate companions to solar-type 



stars both in single and multiple star systems (e.g., Eggenberger et al. 2004). Most 



of these planets have Msin? comparable to or greater than Jupiter, where M is the 



mass of the planet and i is the angle between the orbital plane and a plane perpen- 
dicular to the line of sight to the planetary system. However, as Doppler precision 
increases, lower mass planets are also being discovered and there are now five plan- 
ets with an Msin? less than the mass of Neptune (Extrasolar Planets Encyclopaedia, 
http : I No . ob spm . fr/exoplanetes/ency clo/index .php) . 

Progress is also being made in the development of improved algorithms for the 
analysis of these data. In Gregory (2005a) we described a Bayesian Markov chain Monte 
Carlo (MCMC) approach and demonstrated its effectiveness in the analysis of data from 
HD 73526 (Tinney et al. 2003). This approach has the following advantages: 

a) Good initial guesses of the parameter values are not required. A MCMC is capable 
of efficiently exploring all regions of joint prior parameter space having significant 
probability. There is thus no need to carry out a separate search for the orbital 
period(s). 

b) The method is also capable of fitting a portion of an orbital period, so search periods 
longer than the duration of the data are possible. 

c) The built-in Occam's razor in the Bayesian analysis can save considerable time 
in deciding whether a detection is believable. More complicated models contain 
larger numbers of parameters and thus incur a larger Occam penalty, which is 
automatically incorporated in a Bayesian model selection analysis in a quantitative 
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fashion. The analysis yields the relative probability of each of the models explored. 

d) The analysis yields the full marginal posterior probability density function (PDF) 
for each model parameter, not just the maximum a posterior (MAP) values and a 
Gaussian approximation of their uncertainties. 

e) The inclusion of a novel statistical control system automates the otherwise time 
consuming process of selecting an efficient set of Gaussian parameter proposal 
distribution c's for use in the MCMC algorithm. 

In this paper we will describe a significant improvement to this algorithm and demon- 
strates its use in the re-analysis of data from HD 208487 obtained by Tinney et al. (2005). 
See Ford (2004) for an interesting application of MCMC methods to the interacting plan- 
etary system, GJ 876, in which mutual perturbations are detectable that can be modeled 
with two precessing Keplerian orbits. 

ALGORITHM IMPROVEMENT 

Basically the algorithm was designed as a Bayesian tool for nonlinear model fitting. It is 
well known that nonlinear models can yield multiple solutions (peaks in probability) 
in the prior parameter space of interest. For the Kepler model with sparse data the 
target probability distribution can be very spiky. This is particularly a problem for the 
orbital period parameter, P, which spans the greatest number of decades. In general, the 
sharpness of the peak depends in part on how many periods fit within the duration of 
the data. The previous implementation of the parallel tempering algorithm employed a 
proposal distribution for P which was a Gaussian in the logarithmic of P. This resulted in 
a constant fractional period resolution instead of a fixed absolute resolution, increasing 
the probability of detecting narrow spectral peaks at smaller values of P. However, this 
proved not to be entirely satisfactory because for the HD 73526 data set (Tinney et al. 
2005) one of the three probability peaks (the highest) was not detected in two out of five 
trials (Gregory 2005a). 

The latest implementation employs a proposal distribution in discrete frequency space 
with a frequency interval Av given by Av x data duration = resolution (typically reso- 
lution = 0.007). The resolution corresponds to the fraction of the trial period that the 
last data point will be shifted by a change in search frequency of Av. The width of any 
probability peak is directly related to this quantity so that a proposal distribution of this 
form corresponds more closely to achieving a uniform number of samples per peak. Us- 
ing this search strategy on the HD 73526 data set all three peaks were detected in five 
out of five trials. 



RE-ANALYSIS OF HD 208487 

With the exception of the change mentioned in the previous section, the de- 
tails of the algorithm are as described in Gregory (2005a). Four different models 
(M2j,M2,M\j,&, M\ s ) were compared. Mij& M2j are one and two planet models which 
in addition to the known measurement errors include an empirical estimate of stellar 
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TABLE 1. Parameter prior probability distributions. 


Parameter 


Prior 


Lower bound 


Upper bound 


Pi&P 2 (d) 


Jeffreys 


0.5 


4422 


K x & K 2 (m s- 1 ) 


Modified Jeffreys* 


(K a = 1) 


400 


V (ms- 1 ) 


Uniform improper prior 


all reals 




e\ & e 2 


Uniform 





1 


X1&X2 


Uniform 





1 


0)\ & 0>2 


Uniform 





2k 


s (m s _1 ) 


Modified Jeffreys 


(s a = 5, typical measurement error) 


100 



* Since the prior lower limits for K and s include zero, we used a modified Jeffrey's prior of the form 



p{K\Model,l) = — L- * (1) 

A -+- A fl J n I Aa +Amax j 
\ K a J 

For K <C K a , p(K\Model,I) behaves like a uniform prior and for K ^> K a it behaves like a Jeffreys prior. 
The In ( A: "+5 m " j term in the denominator ensures that the prior is normalized in the interval to K mia . 



jitter which is due in part to flows and inhomogeneities on the stellar surface (e.g., 
Wright 2005). For HD 208487 this is estimated to be somewhere in the range 2 — 4 m 
s _1 . We used a value of 3 m s _1 . In model M\ s we neglect the stellar noise jitter but 
incorporates an additional noise parameter, s, that can allow for any additional noise 
beyond the known measurement uncertainties l . For example, suppose that the star 
actually has two planets, and the model assumes only one is present. In regard to the 
single planet model, the velocity variations induced by the unknown second planet acts 
like an additional unknown noise term. Stellar jitter, if present, will also contribute to 
this extra noise term. In general, nature is more complicated than our model and known 
noise terms. Marginalizing s has the desirable effect of treating anything in the data 
that can't be explained by the model and known measurement errors as noise, leading 
to the most conservative estimates of orbital parameters (see Sections 9.2.3 and 9.2.4 
of Gregory (2005b) for a tutorial demonstration of this point). If there is no extra noise 
then the posterior probability distribution for s will peak at s = 0. 

The parameter prior probability distributions used in this analysis are given in Ta- 
ble 1. For the V parameter, we use a uniform prior which spans all real values and set 
p(V\Model,I) = 1. While this is an 'improper' (non-normalized) prior, the normaliza- 
tion cancels when comparing different models, since it appears exactly once in each 
model. Note, that during the MCMC iterations there is nothing to constrain Pi to be less 



In the absence of detailed knowledge of the sampling distribution for the extra noise, other than that 
it has a finite variance, the maximum entropy principle tells us that a Gaussian distribution would be 
the most conservative choice (i.e., maximally non committal about the information we don't have). We 
will assume the noise variance is finite and adopt a Gaussian distribution with a variance s 2 . Thus, the 
combination of the known errors and extra noise has a Gaussian distribution with variance = of + s 2 , 
where c, is the standard deviation of the known noise for i" 1 data point. 
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TABLE 2. Bayes factors for model comparison. 



Model 


Mean s 


p(D|Model,7) 


Bayes factor 




(m s^ 1 ) 




My 




2.4 ±0.1 x 10" 55 


1.0 


M U 


8 


7.0 ± 0.4 x 10~ 56 


0.3 


M 2 




6.5 ± 0.7 x 10- 54 


27.0 


M 2j 




2.2 ± 0.2 x 10" 54 


9.1 



TABLE 3. Model parameter estimates 



Parameter 


Tinney et al. 


New 

*i n *i 1 a 7 ti t 

diid.1 y »i» 


Parameter 


New 


Pi(d) 


130±1 


129.5l°j 


ft(d) 


998lg 


Ki (ms- 1 ) 


20±2 


16.0+U 


^(ms- 1 ) 


9.8^;9 


<?l 


0.32±0.10 


0.22+;' 4 




0.19±;?| 


©i (deg) 


126 ±40 


126±g 


0)2 (deg) 


157+112 


a x (AU) 


0.49 ±0.04 


0.492t;ooi 


a 2 (AU) 


1 92+ 07 


Mi sin; (Mj) 


0.45 ±0.05 


0.37l;g 


M 2 sin i (My) 


0.46^;?5 


Periastron passage 1 
(JD - 2,450,000) 


1002.8 ±10 


1010^ 


Periastron passage 2 
(JD - 2,450,000) 


4461^ 


RMS residuals (m s -1 ) 


7.2 


4.2 






Reduced % 2 


1.27 


0.83 







than P2- After the MCMC has completed the parameters are redefined so that Pi < P2. 

Table 2 shows the Bayes factors for the 4 models. The Bayes factor is just the ratio of 
the global likelihood of a particular model to model M\j, taken as a reference. Section 
5.3 of Gregory (2005a) describes how the Bayes factors are calculated. Assuming equal 
prior probability for the one and two planet models then Table 2 indicates that a two 
planet model is favored over a one planet model. 

The derived two planet model parameter values are given in Table 3 and compared 
with the one planet parameters given in Tinney et al. (2005). Figures 1 and 2 show 
results for the most favored model M%. The upper panel of Figure 1 shows the raw data, 
the middle panel shows the best two planet fit and the lower panel shows the residuals 
which have an RMS = 4.2 m s _1 . Panel (a) of Figure 2 shows the data, minus the best 
fitting P2 orbit, for two cycles of Pi phase. The best fitting Pi orbit is overlaid. Panel (b) 
shows the data plotted versus P2 phase with the best fitting Pi orbit removed. 
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FIGURE 1. The raw data is shown in panel (a) and the best fitting two planet (Pi = 130 day, P2 = 998 
day) model versus time is shown in (b). Panel (c) shows the residuals. 
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FIGURE 2. Panel (a) shows the data, with the best fitting P2 orbit subtracted, for two cycles of Pi phase 
with the best fitting Pi orbit overlaid. Panel (b) shows the data plotted versus P2 phase with the best fitting 
Pi orbit removed. 



CONCLUSIONS 

In this paper we presented an improved version of the Bayesian parallel tempering 
MCMC algorithm of Gregory (2005a), which is a powerful Kepler periodogram for 
detecting multiple planets in precision radial velocity data. It can easily be modified to 
incorporate precessing Keperian orbits. It is also capable of fitting a portion of an orbital 
period, so early detection of periods longer than the duration of the data are possible. 
From a re-analysis of the HD 208487 data of Tinney et al (2005), we find strong evidence 
for a second planet with a period of 998^ days, an eccentricity of O.^l^g, and an 
Msim = 0.46°^M/. 
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